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ABSTRACT 

This paper presents an efficient method to compute the max- 
imum likelihood (ML) estimation of the parameters of a complex 
2-D sinusoidal, with the complexity order of the FFT. The method 
is based on an accurate barycentric formula for interpolating band- 
limited signals, and on the fact that the ML cost function can be 
viewed as a signal of this type, if the time and frequency variables 
are switched. The method consists in first computing the DFT of the 
data samples, and then locating the maximum of the cost function by 
means of Newton's algorithm. The fact is that the complexity of the 
latter step is small and independent of the data size, since it makes 
use of the barycentric formula for obtaining the values of the cost 
function and its derivatives. Thus, the total complexity order is that 
of the FFT. The method is validated in a numerical example. 

Index Terms — Frequency estimation, parameter estimation, fast 
Fourier transform (FFT), barycentric interpolation, sampling. 

1. INTRODUCTION 

The estimation of the parameters of a complex sinusoidal in either 
one or two dimensions is a pervasive problem in signal processing. 
For this problem, the maximum likelihood (ML) principle provides 
an optimal estimator, in the sense that it achieves the Cramer-Rao 
bound at relatively low signal-to-noise ratios. Yet in practice this es- 
timator is deemed too complex, due to the associated maximization 
problem. In rough terms, the situation is that, though the ML cost 
function can be regularly sampled using the FFT efficiently, the lo- 
calization of its maximum requires some search procedure, and here 
is where the high computational burden seems unavoidable. This ob- 
servation was first stated by B. G. Quinn in [ 1 1, where it was shown 
that a Gauss-Newton iteration may fail to find the global maximum 
of the cost function, if the initial iteration is taken from the DFT 
of the data samples (without zero padding). This has led several 
researches to give up the ML approach, and look for sub-optimal es- 
timators with lower computational burden, J2]|3]- However, in some 
references [4. 5 6 1 there has been an attempt to recover the initial 
ML approach, based on two arguments. The first is that it is feasi- 
ble to approximately locate the maximum of the ML cost function, 
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simply by looking for the DFT sample with largest module. Be- 
sides, the accuracy of this coarse localization can be improved by 
zero padding. And the second is that it is possible to interpolate the 
cost function close to this location from the nearby DFT samples 
with some accuracy, since it is a smooth function. These arguments, 
if properly exploited, make it possible to improve the performance 
significantly, with a complexity similar to that of the FFT. The pur- 
pose of this paper is to go one step further in this direction, and 
show that it is feasible to perform this interpolation with high accu- 
racy from a small number of DFT samples, so that the actual ML 
estimate can be obtained with the complexity of a single FFT. The 
key lies in viewing the DFT as a band-limited signal in the frequency 
variable. 

The basic concept in this paper is first presented in the next sec- 
tion in the context of time-domain interpolation. It is a simple an 
efficient technique to interpolate a band-limited signal from a small 
number of samples with high accuracy. Then, the estimation prob- 
lem in one dimension is introduced in Sec. [3] where it is shown that 
the method in the next section is actually the key for an efficient solu- 
tion, if the independent variable is properly interpreted. Afterward, 
Sec. |4]presents the extension of the method in Sec. [3]to the two- 
dimensional case, and finally Sec. [5]contains a numerical example. 

2. HIGH ACCURACY INTERPOLATION OF A 
BAND-LIMITED SIGNAL AND ITS DERIVATIVES 

Given a band-limited signal s(t) with two-sided bandwidth B, a 
usual task is to interpolate its value from a finite set of samples sur- 
rounding t. If the sampling period is T with BT < 1 (Nyquist 
condition), and 2P + 1 samples are taken symmetrically around t, 
then this task consists in finding a set of coefficients a p (u) such that 
the formula 

p 

s(t) « Y, s((n - p)T)a P (u) (1) 
P =-p 

is accurate, where t — nT + u is the modulo-T decomposition of t, 



[t/T + 1/2J, 



t - nT. 



(2) 



For fixed it, the a p (u) can be obtained numerically using filter op- 
timization techniques |7). Yet this approach is cumbersome, if not 
only s(t) but also its derivatives must be interpolated for varying 
values of u efficiently. Recently, a so-called barycentric interpola- 
tor was derived in [8] that solves this problem satisfactorily. This 
interpolator takes the form 



s(nT + u)^ ^ 
P =-p 



s((n — p)T)w p 
u-pT _ 



E • 

P =-p 



pT 



(3) 



where w p is a set of constants which are samples of a fixed function 
w(t), [w p = w(pT)]. In [8|, this function is given by the formula 

M = r(t/r + p + i)r(-t/r + p+i)g(t) 

L'(t) ' (> 

where T(-) is the Gamma function, g(f ) is the pulse 



8inc((l-BT)V(t/r)a-(P+l)2) 
sinc(j(l-BT)(P+l)) 

and L(£) is the polynomial 

p 

L(t)= n 
P =-p 



(5) 



(6) 



(See Q0 for further details.) The fact is that the error of ([T} decreases 
exponentially with trend e~ 7r( - 1 ~ BJ " >p . In practice this means that a 
small P is enough to obtain high accuracy. Besides, as shown in [8|, 
Eq. l[T) can be differentiated with low complexity, so as to interpolate 
the differentials of s(i) of any order, and ((3} can be evaluated using 
only one division. This interpolator will be the fundamental tool in 
the next section. 

3. ML FREQUENCY ESTIMATION: 1-D CASE 

Consider a signal z(t) consisting of an undamped exponential with 
complex amplitude a and frequency f , contaminated by complex 
white noise n(t) of zero mean and variance a 2 , 



z(t) 



J** to* 



+ n(t). 



(7) 



For simplicity, f is assumed to lie in [—1/2, l/2[. Next, assume that 
M samples ofz(i) are taken at instants mi, mi + 1, . . . ,mi+M— 1 
with mi = — [M/2]. The maximum likelihood estimation of /„ 
from these samples is the argument that maximizes the cost function 



L(/) = |c(/)| 2 , 
where c(/) is the correlation 



mi+M-l 



c(/) 



E 



-j2wj'm 



(8) 



(9) 



t in l[T). Also, L(/) is band-limited with bandwidth 4mi. So, the 
method of sampling the correlation using the FFT, and then looking 
for its maximum module can be reformulated taking into account 
these information. First, since L(/) has two-sided bandwidth 4mi, 
it is necessary to sample this function with a spacing smaller than its 
Nyquist period l/(4mi), in order to coarsely locate its maximum. 
This implies that a factor-two zero padding is enough to ensure the 
localization of the maximum of L(/). And second, since c(/) has 
bandwidth 2m i, it can be interpolated using the bary centric formula 
in Q} with 2mi in place of B, and Af in place of T, i.e, 



=(/)* E 



P =-p 



c((k~p)Af)w p 
7 ~ PA/ 



E 



j-pAf 



(11) 



where / = kAf + 7 is the modulo-A/ decomposition of / de- 
fined as in l(2j. If £(/) denotes the approximation to c(/) in i ll It . 
then one may replace the problem of maximizing L(/) with that of 
maximizing 



Hf) = i5(/)r 



(12) 



Besides, the maximum of this function is close to the abscissa of the 
largest FFT sample in dlOt , k Af. 

Since P is small and there is a coarse estimate of the maximum 
abscissa, the maximization of j 1 2b is a low-complexity problem that 
can be solved using standard numerical methods. Given that the 
differentials of the barycentric interpolator can be easily computed 
(8] Sec. IV], a suitable numerical method is Newton's algorithm, in 
which the rth iteration f r is refined using 



fr 



/r + L'(/r)/L"(/ r ) 



where 



£'(/) = ^Re{c'(/)c(/)*}, 



£"(/)= 



A(Re{c"(/)c(/)*}+|c'(/)| 2 ) 



(13) 



(14) 



This process can be initiated with /1 = k Af, and a small number 
of iterations (3 to 5) is enough to obtain f . 

The computational burden of this method is given by that of the 
FFT, i.e, it is 0{M log M), and it yields the actual ML estimate of 

U 



Since c(/) is the DFT of the samples z(m), the maximum of L(/) 
can be approximately located by selecting a frequency spacing Af 
such that 1 / Af is a power of two, and then computing the samples 



mi -j-M- 

;(fcA/) = ]T 



z(m)e 



-j2ir fcmA/ 



(10) 



by means of a radix-2 FFT algorithm. The cost of this operation is 
only 0(Af logM). This way, it is possible to obtain a frequency 
k Af that lies close to the true abscissa of the maximum of L(/). 
However, if f a is this abscissa, the approximation of f„ using k A f 
is very inaccurate. In this situation, the accuracy can be improved 
by reducing Af (over-sampling), but then the computational burden 
becomes high. 

The barycentric interpolator in the previous section provides an 
efficient method to obtain f a The key idea is that the correlation 
c(/) in ® is a band-limited signal in the f variable of bandwidth 
2mi, i.e, the variable / in ((9} plays the same role as the variable 



4. ML FREQUENCY ESTIMATION: 2-D CASE 

The method in the previous section can be extended to two- 
dimensional signals with non-essential changes. For this, it is 
only necessary to consider two variables, t\ and £2, and repeat the 
same interpolation procedure. The initial model, equivalent to ((7), is 



z(ii,*2) 



pj'2ir(/l,otl-|-/2, t2) 



+ n(t u t 2 ). 



(15) 



Next, this signal is sampled at the integer pairs (m, n) for mi < 
m < mi + M and n\ < n < m + N, with mi = - [M/2] , 
m = — [iV/2] . The 2-D equivalent of the cost function in {8]l is 



L(A,/ 2 ) = |c(/i,/ a )| 
where c(/i, /a) is the correlation 

m-i+M— 1 ni+JV-1 



c(/l,/2 



E E z ( m . n ) 



-j2-n-(f 1 m+f 2 n) 



(16) 



(17) 



This function can be sampled with spacings A/i, A/ 2 using a 
radix-2 2-D FFT. These samples provide a frequency pair (fci, A/i , 
fe.oA/j) that lies close to the maximum of L(/i,/2). Then, it 
is possible to set up an interpolation formula like ( lilt but in two 
dimensions, 



Pi Pi 

c(A,/ a )» E E 

Pl=-Pl P2 = -f2 

c((fci -pi)A/i, (fe -p 2 )A/ 2 )wi, pl 'u;2,p 2 
(71 — piA/i)(7i -piA/i) 

L ^ 71-P1A/1 ^ 72 -p 2 A/ 2 J ' 

Pl=--Pl vi = -Pi 



(18) 



where /1 = feiA/i + 71 and / 2 = fc 2 A/ 2 + 7 2 are the modulo 
decompositions, and Wi lP1 and «; 2jP2 are the barycentric weights 
corresponding to bandwidths 2mi and 2ni, and truncation indices 
Pi and P 2 , respectively. If c(/i, / 2 ) denotes the formula in ([18}, the 
problem of maximizing L(/i, / 2 ) can be substituted by the problem 
of maximizing 

t(/i,/ 2 )s|c(/i,/ a )| 2 . (19) 
Finally, the Newton iteration for this problem is 



C/V+1,/2,,+1) = (/!,,., / 2 , r ) + WL _1 VL, 



(20) 



where VL and T-CL are the gradient and Hessian of L respectively, 
evaluated at (fi, r , / 2 ,r). These functionals are 



VL : 



MN 



Re{Vcc*} 



and 



MN 



Re{Ucc* + VcVc H }, 



(21) 



(22) 



where Vc and He are the corresponding gradient and Hessian of 5. 

Since the evaluation of c(/i, / 2 ) has a small cost which is inde- 
pendent of either M or N, the complexity is given by the 2-D FFT, 
i.e, it is 0(MNlog(MN)), Note that sub-optimal methods like 
that in [3| have complexity 0(MN(M + N)), that is, the method 
in this section yields the ML estimate and, besides, has a smaller 
complexity order. 

5. NUMERICAL EXAMPLES 

Let 4>(f, t) denote the value delivered by the barycentric formula 
in l[3}, when the input signal is the undamped exponential s(i) = 
e j27r/t § mce me f unc tions to interpolate in this paper are sums of 
exponentials like e- ?2,r ^', a simple way to assess the interpolation 
accuracy is to evaluate the maximum module of e? 2 ^f u — (j>(f, it), 
for |u| < T/2 and / varying in [~B/2, B/2], i.e, to assess the 
spectrum function 



E(/) 



max 

u\<T/2 



j2irfu 



(23) 



FigQ]shows E(/) for BT = 0.25 and several truncation indices P. 
Note that any accuracy can be achieved uniformly in [—B/2, B/2], 
by slightly increasing P. 

In order to test the interpolation error in a specific example, a 
2-D undamped exponential with M — 500, iV = 651 and frequen- 
cies / ,i = 0.234452 and / , 2 = -0.143254 was generated. Then, 
complex white noise was added, so that the signal-to-noise ratio was 




0.02 0.04 0.06 0.08 0.1 0.12 
Normalized frequency f/T 

Fig. 1. Error spectrum E(/) of the barycentric interpolator in ((3] for 
several truncation indices P. 
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Fig. 2. Truncation index P versus the log 10 of the maximum inter- 
polation error (log 10 e) for the cost function L(/i, / 2 ). 



S/N = 5 dB. Fig.[2]shows the interpolation error for the cost func- 
tion L(/i, / 2 ), defined by 



max /i,J 2 m/1,/2) -L(/i,/ 2 )| 
max fuf2 L(/i,/ 2 ) 



(24) 



for varying truncation index P. Again, it is clear that any accuracy 
can be achieved by slightly increasing P. 

Next, two estimators were compared in this example. To de- 
scribe them, let Z denote the data matrix obtained by sampling d!5t 
as described in that section, 



z((mi + m - l)Ti, (m + n - l)T a ), 



(25) 



and let Ui, vi denote its first left and right singular vectors, respec- 
tively. The first method consisted in computing the 1-D ML estima- 
tor from ui and vi so as to obtain the respective estimations of f 0i i 
and f , 2 , 



-j2ir(m 1 —l + m)f 



/ ,i = argmax j ^ [ui] m e 

m — 1 
N 

/V V i — j2ir(ni — 1+n 
0,2 = argmax | ) Jvij^e 



(26) 
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Fig. 3. Root-Mean-Square (RMS) error of the subspace and ML 
estimators, together with the Cramer-Rao bound. 
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Fig. 4. Average number of iterations required by the Newton 
method. 



The actual estimations f ,i and / 0i 2 were computed from ui and 
Vi by applying the 1-D method in Sec. [5] to each of them. This 
method had complexity 0(MN(M + N)), since it involved the 
computation of the singular vectors ui and Vi, and is equivalent 
to the method in |3|. The second method computed the actual ML 
estimator from Z using the method in Sec. [4] and its complexity was 
just 0(MNlog(MN)). For the values of M and N given above, 
the over-sampling factors were 2.048 and 1.523, respectively, i.e, 
the FFT had size 1024 in both dimensions. 

Fig. [3]shows the root-mean-square (RMS) error for the first es- 
timator, termed subspace estimator, for the second estimator (inter- 
polated ML estimator), and the Cramer-Rao bound. The threshold 
performance of the second is clearly superior. 

Finally, Fig. [4] shows the average number of iterations that re- 
quired the interpolated ML estimator, which was roughly equal to 
three. 



high accuracy from a small number of samples, if the sampling fre- 
quency is somewhat higher than the Nyquist frequency. Besides, a 
specific barycentric formula is proposed to perform this kind of in- 
terpolation. And second, it is shown that the ML cost function for the 
estimation of a complex 2-D (and 1-D) sinusoidal can be viewed as a 
band-limited signal, if the time and frequency variables are switched. 
Finally, these two results are combined in a method that is able to de- 
liver the ML estimate with the complexity order of the FFT, which 
is based on Newton's algorithm. 
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6. CONCLUSIONS 

A method has been presented that allows one to compute the maxi- 
mum likelihood (ML) estimation of a complex 2-D sinusoidal, with 
the complexity of the fast Fourier transform (FFT). First, it is re- 
called in the paper that a band-limited signal can be interpolated with 



